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Abstract 

£H ■ We consider the spreading of a thin two-dimensional droplet on a planar substrate as 

a prototype system to compare the contemporary model for contact line motion based 
on interface formation of Shikhmurzaev [Int. J. Multiphas. Flow 19, 589 (1993)], to the 
more commonly used continuum fluid dynamical equations augmented with the Navier- 
slip condition. Considering quasistatic droplet evolution and using the method of matched 
asymptotics, we find that the evolution of the droplet radius using the interface formation 
model reduces to an equivalent expression for a slip model, where the prescribed micro- 
scopic dynamic contact angle has a velocity dependent correction to its static value. This 
^>~ > ' result is found for both the original interface formation model formulation and for a more 

tIh ■ recent version, where mass transfer from bulk to surface layers is accounted for through 

the boundary conditions. Various features of the model, such as the pressure behaviour 
and rolling motion at the contact line, and their relevance, are also considered in the 
prototype system we adopt. 

1 Introduction 

A contact line is formed wherever two immiscible fluids, such as a liquid and a gas (often air) 
meet a solid surface. The wetting of a solid by a liquid occurs via a moving contact line in 
the gas/liquid/solid system, and is a crucial part of many natural and technological processes 
(see e.g. the review articles of de Gennes, 1 Blake, 2 and Bonn et al.—). A model developed 
by Shikhmurzaev, 4 building on the earlier work of Bedeaux et al.— on non-equilibrium ther- 
modynamics for immiscible fluids, is one of many put forward in the literature to understand 
this apparently simple situation. Since its inception the model has been used to describe a 
more general class of problems where interface formation occurs,— and as such, is termed the 
"interface formation model". 

The application of the usual hydrodynamic equations and fluid-mechanical modelling as- 
sumptions to the moving contact line leads to the nonexistence of a solution.— This issue was 
identified by Moffatt, 8 and Huh and Scriven— for a planar interface leading to the well-known 
problem of a nonintegrable shear-stress singularity. To alleviate the singularity, a variety of 
models and associated contributing phenomena have been proposed. These include allowing 
slip to occur at the solid surface, a rather popular model; accounting for the existence of a pre- 
cursor film ahead of the contact line,— a popular model also; incorporating rheological effects 
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such as shear-thinning;— incorporating intermolecular forces and considering the interface to 
be diffuse;— or including evaporative fluxes 1 ^ (for a review of some of these approaches in 
the context of the interface formation model see §3.4 of the monograph of Shikhmurzaev 6 ). 
Savva and Kalliadasis^ have recently demonstrated an equivalence between a selection of 
some of the most popular models, i.e. slip and precursor film models, in spreading situations 
(building on the work of Pismen and Eggers-i»), so here we primarily assess the interface 
formation model in comparison to the popular slip model mentioned earlier, the Navier-slip 
condition 1 — detailed in Sec. 11.11 

A feature, reported to be significant, 16 of the interface formation model over the usual 
hydrodynamic equations augmented with the Navier-slip boundary condition is that the mi- 
croscopic dynamic contact angle (as opposed to the apparent contact angle at observable 
distances away from the contact line) is determined as part of the solution. That no em- 
pirical relation between contact line velocity and contact angle is required, and additionally 
that the influence of the flow field on the contact angle is incorporated, makes the interface 
formation model worthy of analysis. It should be noted, however, that this microscopic con- 
tact angle is determined through contact line conditions arising from a mechanical balance 
local to the contact line for the interface formation model, where similar reasoning leads to 
the imposition of the static Young's contact angle for the Navier-slip model at all times. 
Empirical relations may also be applied for the microscopic contact angle variation with slip 
models,— but it is the fact that this is not required for the interface formation model that is 
of interest — especially as empirical relations will be based on experiments at observable dis- 
tances from the contact line, but then enforced for the contact angle at the microscale. It is 
also worth mentioning here that precursor film models, where a very thin film covers the solid 
surface ahead of the macroscopic liquid/gas interface through the inclusion of a disjoining 
pressure (due to long-range intermolecular interactions), remove the contact line singularity 
as there is no longer a sharp contact line, but an apparent one. In these models this apparent 
contact angle is related to the model parameters, such as the Hamaker constant and film 
thickness, but no actual microscopic contact angle is formed. For further details of these and 
related aspects of the interface formation model, we refer the reader to Shikhmurzaev's earlier 
works 4,16 and monograph. 6 

However, as one might expect, the model has not remained free of controversy. There have 
been occasions where other authors have questioned its experimental validity (see Lindner- 
Silwester and Schneider, 1 -- noting the response of Shikhmurzaev— and their reply—), the 
assumptions underlying its formulation (see Eggers and Evans^i and the response^), and 
the possibility of the requirement of an additional boundary condition at the contact line for 
consistency (suggested by Bedeaux— and Billingham 24 ). In a recent issue of the European 
Physical Journal - Special Topics, it is clear that the debate around the interface forma- 
tion model remains very active, with strong opinions from a variety of authors from diverse 
scientific backgrounds published there.— ~— 

The two papers by Billingham,— >2I mentioned above, for the steady motion of the contact 
line in a wedge and for a thin-film moving under gravity, give arguably the most careful 
analytical work with the interface formation model. We will consider the model in a similar 
framework but for a different problem, that of quasistatic droplet spreading. This prototype 
system will also help to provide insight into the interface formation model, as we intend to 
scrutinise the model in light of the above opinions and debate. 

Another feature of the model, much lauded by its advocates, is that in a liquid/gas/solid 
system, the liquid is predicted to display rolling motion. This flow pattern has been experi- 
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mentally observed by a variety of authors (see for instance Dussan V. and Davis^ and Chen 
et a/., 39,40 the latter with an inner limit at ~ 30 |J.m from the contact line), but observations 
are on the micro- rather than the nano-scale, where significant differences between models 
will appear. Rolling motion is built into the interface formation model through its formula- 
tion, whereas in the "classical model" the contact line instead slips over the solid. Whilst 
this point will be evident from our analysis, and it is mentioned as a benefit of the interface 
formation model over others by Shikhmurzaev — we reiterate that experimental observations 
of flow patterns at the microscale are in a region where slip models will also show similar flow 
patterns — albeit with the possibility of the stagnation point causing a slowing of the flow in 
the vicinity of the contact line for slip models. We also highlight the experimental work of 
Savelski et al.,— where it is suggested that rolling motion may instead occur in the gas phase 
with a splitting streamline present in the liquid (for particular values of the apparent dynamic 
contact angle and the viscosity ratio between the two fluids). These results are questioned by 
Shikhmurzaev— ^ suggesting that the Reynolds number in such experiments must be based 
on the observable distance from the contact line and that the results do not invalidate the 
rolling pattern near the contact line where the Reynolds number tends to zero. 

In earlier works of Shikhmurzaev (e.g.— ^) and in the appraisals of Billingham, 24 ' 37 mass 
exchange between the interfaces and the bulk was taken into account in the surface phase 
mass balance equations but neglected in the boundary conditions for the normal compo- 
nent of the bulk velocity. The neglected terms were found to be critical to the convergent 
flow of a Newtonian fluid near a free boundary,— and are included in the version given in 
Shikhmurzaev's monograph. 6 It is claimed that the earlier works considering dynamic wetting 
at small capillary numbers do not need to be revisited in light of the additional mass transfer 
contributions.— «2 Interestingly however, when considering the region close to the contact line 
we find in Sec. 13.2.21 that the low-velocity region predicted by slip models and preventing 
nanoscale rolling motion also occurs for the original formulation of the interface formation 
model. The additional terms in the modern formulation overcome this low-velocity region en- 
abling rolling motion to take place, seen in Sec.[H By considering the influence of the imposed 
interfacial boundary conditions, it may be concluded that a stagnation point at the contact 
line will occur for any model where mass transfer is not allowed. Diffuse interface methods, 
through the creation of an interfacial region of density variation, include mass transfer and 
should thus allow nanoscale rolling motion. A connection between diffuse interface methods 
and the interface formation model has been briefly suggested previously, 27 ' 43 but a detailed 
comparison would be of interest (although it lies beyond the scope of the analysis presented 
here) . 

Whilst there have been a number of notable uses of the interface formation model, 42l44 ~ — 
there is still much room for further work. That there is not a much greater use of the model in 
the literature may in part be due to its comparative complexity, although it is also likely that 
the controversy discussed above alongside questions over its physical basis have dissuaded 
others to attempt detailed independent analyses of the model. Nevertheless, whilst we have 
attempted to highlight the ongoing debate, we believe that the model does have a number of 
advantages, and that further understanding of both the model and the questions raised about 
its formulation could aid in the long-standing issue of the moving contact line problem. 

It is the intention of this work to consider the interface formation model in the simple 
yet illustrative situation of two-dimensional droplet spreading on a flat horizontal substrate 
utilising the long- wave approximation appropriate for thin droplets. The spreading of a 
droplet has been widely considered as a paradigmatic system to evaluate dynamic contact 
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line motion by many authors. In particular, Hocking considered the competition of gravity 
and capillarity with Navier-slip,— and Haley and Miksis extended this for other similar slip 
models.— The interface formation model has also been considered for droplet spreading— 
with an extension for microdrops,— but not within the long-wave approximation as is the 
case here. More specifically, by using this approximation for the problem of droplet spreading 
on a planar horizontal substrate, as in the framework of Hocking,— we obtain a substantially 
reduced setting which is amenable to mathematical and numerical scrutiny - thus allowing 
for direct comparisons between Navier-slip and interface formation models to be made with 
relative ease. 

This work is outlined as follows. In the remainder of this section (in ll.lHl.2p we reprise 
the classical model and provide a summary of the interface formation model with some obser- 
vations on the main points of its formulation. Then in Sec. 12. II we detail the model equations 
and their forms in our droplet regime, followed by the nondimensionalization and long-wave 
assumption in Sec. 12.21 A matched asymptotic analysis is performed in Sec. [3] for the original 
interface formation model, with changes to the analysis for the modern formulation in Sec. [H 
Our asymptotic matching procedure is verified through considering numerical solutions to the 
full problems, details in Sec. [5j with conclusions in Sec. [6j 



1.1 Classical hydrodynamic model 

For the classical model, the conditions on a flat solid boundary with unit normal and tangent 
vectors rig, t$ are, no-slip, u • = 0, and no-penetration, u • rig = 0, where u is the bulk 
fluid velocity. As detailed earlier, enforcing the no-slip condition leads to the nonexistence of 
a solution to the moving contact line problem. To alleviate this difficulty the classical model 
invokes the popular Navier-slip condition 

u-tg= ns-er-ts, (1) 

A* 

where \i is the fluid viscosity, /3ns is the coefficient of slip, and a = — PI + T is the total stress 
tensor, where P is the pressure, I is the identity tensor, and T = 2^iD = /i [(Vu) + (Vu) T ] 
is the extra-stress tensor. On a free surface between a liquid and a gas with unit normal and 
tangent vectors ng, t G , we have 

^ = |^ + u-V/ = 0, n G -T-t G = -V s a NS -t G: P G -P + n G -T-n G = ajv 5 V s -n G , 

(2) 

which are the kinematic condition, the tangential stress equation, and the normal stress 
equation, where ctns is the surface tension on the liquid-gas interface f(r,9) = 0, P G is the 
pressure of the gas, and T is time. The tangential stress equation gives rise to the Marangoni 
effect when a surface tension gradient is present, caused by gradients in temperature or 
variation in chemical composition at the interface. V s ctns is the surface gradient of <jns- F° r 
scalar field / and vector field f, we define 

V s / = (I-n®n).V/, V s .f=tr((I-n®n) • Vf) , (3) 

where on the right hand side, V is the usual gradient operator, and n (8) n is the dyadic 
product of the normal vector n with itself.— 

Finally, at the contact line there is no way of determining the microscopic dynamic contact 
angle 9d from the classical model. As such it must be prescribed to provide a boundary 
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condition at the contact line, either by assuming that it is always equal to the static contact 
angle, 9^ = 9 S , or by imposing a relationship on the microscopic dynamic contact angle 
in terms of the contact line velocity. The interface formation model by comparison has the 
microscopic dynamic contact angle determined as part of the solution via the Young equation, 
as mentioned earlier, so that no empirical relation is required. 



1.2 The interface formation model 

There are two interrelated mechanisms key to the formulation of the interface formation 
model: 

1. Surface layers: In a region near to an interface between a liquid and another substance, 
liquid molecules are not fully surrounded by other liquid molecules as they would be 
in the bulk due to the proximity of the other material. As such they lose a number of 
cohesive interactions causing an uneven force on the liquid molecules and giving rise to 
a thin surface layer where the liquid density differs from the bulk (Biilingham 2 - quotes 
a typical thickness of 10 -10 m). This variation in density viewed on a macroscopic scale, 
appears as the surface tension of the liquid, see Fig. [TJ This can also be thought of as 
liquid molecules near the surface being in an unfavourable energy state and the surface 
tension being a direct measure of this energy shortfall per unit area.— 

Consider a surface layer with surface velocity v s . The conditions on either side of the 
layer may differ — for instance the surface layer between liquid and solid may exhibit no- 
slip on the solid- facing side (in accordance with the classical model), but on the liquid- 
facing side may slip, see Fig. [2]Ja). The continuum mechanical limit of these surface 
layers being of zero thickness is then taken, and thus the surface properties influence 
the bulk simply by modifying the boundary conditions applied at each interface. We 
summarise the implication of surface layers as follows: 

• Density in surface layer differs from that of the bulk due to presence of other 
material, 

• Density variation gives rise to surface tension, 

• Allows different conditions to be applied on each side of the surface layer — no- 
slip is then satisfied for the fluid particles actually next to the solid, with a slip 
condition appearing as the boundary condition due to averaging throughout the 
surface layer. 

The surface layers have the properties of surface velocity v s , surface tension a s , surface 
density p s , equilibrium surface density pi and a surface density associated with zero 
surface tension p^. These will be differentiated between the solid- liquid and the gas- 
liquid interfaces with subscripts S and G respectively. 

2. Surface tension relaxation: Consider a static contact line as in Fig. [2](b), where the 
surface tensions associated with the interfaces between fluids 1 and 2 and the solid 
surface are au, ais and o^s- In this situation, the horizontal force balance is given by 
the Young equation: 

(712 COS 6 S + Ois = CT2S- (4) 

[It is noteworthy, that beside the Young equation, which reflects the mechanical equi- 
librium in a direction parallel to the wall at the contact line, there exists a relation in 
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the normal directio n 55 i 56 but its implications are beyond the scope of this study.] For 
a moving contact line also holds as any momentum fluxes due to the motion of the 
interfaces are negligible (see §4.3.5.2 of Ref. @), but in this case with dynamic values of 
the surfaces tensions. When in motion, the microscopic dynamic contact angle varies 
from its static value causing either 012, 0"is, &2S ° r a combination to also deviate from 
their static values. Discussion about the Young equation including some recent disputes 
is also given in §2.4.4 of the monograph of Shikhmurzaev.— 

Material flux through the contact line is caused by the rolling motion of the fluid. 
As the properties of the liquid-gas and liquid-solid interfaces will not be the same in 
general, then as a fluid particle is transferred from the gas interface to become part of 
the solid interface, some reorganisation of the molecules will be required. In particular, 
a fluid particle is associated with the surface tension a\2 when on the gas interface, but 
after moving through the contact line must then be associated with the surface tension 
ais- The crucial assumption here is that this relaxation of the surface tension to the 
equilibrium existing far from the contact line happens in a finite time, r, rather than 
instantaneously,— 1 ^ this assumption being vital for the interface formation model as an 
instantaneous relaxation will reduce the model to that of classical Navier-slip. 
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Figure 1: The density variation close to an interface, and the resulting surface tension. Fluid 
particle A away from the interface is surrounded by other fluid particles, whereas the fluid 
particle B near the interface feels the presence of fewer fluid particles and is thus in an 
unfavourable state. It is this which gives rise to the surface tension on continuum mechanical 
scales. 
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Figure 2: (a) Sketch of the surface layer with velocity at the solid surface satisfying no-slip. 
This surface layer is taken with zero thickness in the continuum approximation so that the 
bulk velocity u undergoes apparent slip at the boundary, (b) A Young's force diagram showing 
the forces acting at a static contact line between two fluids and a solid substrate, with the 
three phases meeting at static contact angle 6$ (see Billingham— ) . 

It is important to point out that in the following we have the unit normal at both solid 
and gas interfaces 115 and pointing into the liquid as defined by Shikhmurzaev,— where 
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we also refer the reader for further information about the interface formation model and for 
details of its derivation. 



The interface formation model equations We now consider the equations of the inter- 
face formation model. First consider the flow in the bulk liquid, which is governed by the 
usual Navier-Stokes equations for a viscous incompressible Newtonian fluid 

V-u = 0, p(J^ + u-Vu) =-VP + /iV 2 u, (5) 

where p is the bulk density, other variables are as for the classical model and gravitational 
effects are assumed to be negligible. We now consider the solid-liquid interface: 

v% • n s = 0, (6) 

t«tf+2 V » <7 5 = J 8u ||» p(u-v|)-n s = Xm Ps ~ PSe , (7a,b) 

^ + V s • (p|v|) = - Ps ~ PSe , v s 5|| = i U| | + aV 8 a%, a s s = 7 (^ 0)5 - p%), 

(8 a,b,c) 

where U|| = (I — ng <8> ns) • u and = (I — n$ ® ng) • v s s are the tangential projections 
of the bulk and surface velocities at the solid surface, and t xy = 2pns • D • (I — ns <8> ns) 
is the shear stress exerted by the bulk fluid on the interface. Selecting \ m = or \ m = 1 
determines whether mass transfer between the surface layer and bulk is included throughout 
the surface layer equations, giving either the original or modern formulation of the interface 
formation model. These equations are for a static solid substrate so there is no solid velocity 
contribution, in comparison to the equations in the monograph of Shikhmurzaev.— We choose 
to work in a frame of reference of a static observer (such as also taken by Billingham— and 
Shikhmurzaev—) . 

The equations above have been grouped to broadly indicate where they apply, if imagining 
the surface with a thickness for illustrative purposes. Equation © applies at the solid-facing 
side of the surface layer, the condition describing the usual no-penetration of the liquid. 
Equations (7) apply at the liquid-facing side, being a generalised Navier-slip condition and 
a condition on the normal component of the bulk velocity (with mass transfer between bulk 
and surface layer included or neglected through the choice of Xm)- Equations (8) apply 
throughout the surface layer, giving the mass balance, a surface velocity relation (being a 
jump in tangential velocities due to the surface tension gradient), and a relationship between 
surface tension and surface density — a linear relationship is assumed, although Shikhmurzaev 
gives suggestions as to how this may be generalised if required. - If there was no surface tension 
gradient, the surface velocity relation would suggest the surface velocity is the average of the 
bulk and solid velocities, whereas with surface tension gradient, the effect is similar to Darcy's 
law for the pressure drop in a porous pipe.— 

We note that a, (3 and 7 are parameters of the surface (typical values quoted in Billing- 
ham—, see also Table 1), with a corresponding to the response of the interface to surface 
tension gradients, f3 to the coefficient of sliding friction and 7 to the inverse surface layer 
compressibility. The partial derivative of a surface variable, such as in (8a), needs to be 
considered with care, but any differences between definitions are negligible in the long-wave 
approximation.— 
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Next we consider the gas-liquid interface: 

f£ = i+v & .v/=o, ^ m 

P G -P + n G -T-n G = a s G V s -n G , t xy = -V s a s G , p(u - v^) • n G = Xm Pc ~ P ° e , 

T 

(10 a,b,c) 

9 -§ + V s • (pfcvfc) = -^L^k, y^Cv^i - u,|) = V s a G , a s G = 7 (^ 0)G - p G ), 

(11 a,b,c) 

where similarly U|| = (I — ® n G ) • u and = (I — <8) n G ) ' v g are the tangential 
projections of the bulk and surface velocities, and t xy = 2^n G • D • (I — <S> n G ) is the shear 
stress. Equation © gives the usual kinematic condition, applying on the gas-facing side of the 
surface, the free surface being associated with the surface velocity v^. Equations (10) apply 
at the liquid-facing side and are the normal and tangential stress conditions (where viscous gas 
stresses are assumed negligible, as in many other works with both interface formation model— 
and classical model&i^^i) , and a condition on the normal component of the bulk velocity. 
Unlike the classical formulation the gradient of surface tension plays a role here purely based 
on the flow of the fluid — with no variation in temperature or chemical composition required 
(Marangoni effect). Equations (11) finally are the analogues of (8) above, with the unusual 
factor in (lib) readily derived from general (rather than specific to solid- liquid or liquid- 
gas) surface equations in the form of (7a), (8b) and (10b), see Ref. for details of the 
derivation. It should be noted here that a, [3, 7 and r have been used for both surface layers 
and for simplicity are assumed to be equal for both (consistent with other works with this 
model^^). 

Finally, we require conditions at the contact line. It is postulated in Refs. 0@ that there 
are two conditions at the contact line. 

• Conservation of mass through the contact line: 

p G (v G - V CL ) ■ t G + p s s (v s s - U CL ) •t s + Q = 0, (12) 

where we have included the contribution from the contact line velocity, TJql since (as 
mentioned earlier) we are in the frame of reference of a static observer. Equation (|12|) 
states that the mass flowing into the contact line from the gas-liquid surface equals the 
amount flowing out into the solid-liquid surface, with the possibility of mass transfer to 
the bulk with non-zero Q. It is assumed in most previous works that Q = 0, and we 
make the same assumption here. 

• A force balance at the contact line, from 

a s G cos d + 0% = a s Ge cos 6 8 + a% & (13) 

where the assumption is made that the surface tension associated with the solid-gas 
interface, a G , does not vary from its static value as the contact line moves. 37149 
Shikhmurzaev comments that a G is negligible for all contact angles except if close to 
180°, leading to the same condition above. 6 As such there appears to be unanimous 
agreement on this boundary condition as we use it. 
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A third condition is found necessary by Billingham 24 and derived by Bedeaux— through 
consideration of the entropy production rate at the contact line, which relies on the ability 
to treat the contact line as a separate thermodynamic system — with associated properties 
such as line tension. Shikhmurzaev, in §4.3.5.3 of his monograph— puts forward his point of 
view regarding this treatment, suggesting that a line phase cannot exist whilst remaining in 
the continuum approximation, referencing other authors such as de Gennes.— There remains, 
however, much debate and disagreement over this point and if one follows the recent work of 
Bedeaux~, then the treatment of the line phase considering conservation of energy arguments 
leads to the condition 



p s s (v s s - U C l) • t s = U (p s s - p s Se ) + U G0 (p s G - p s Ge ), (14) 

where Uq and Uqq are constants with dimensions of velocity and assumed to be equal for 
simplicity (as assumed by Billingham, with (|14p given as found there). This condition 
states that the flux through the contact line is driven by the deviation of the surface layer 
densities from equilibrium. 

As mentioned earlier, a full derivation of the interface formation model equations can be 
found in Ref. @- We highlight that equations ([6]), ([9]) and ([12]) arise through conservation 
of mass, and the normal and tangential stress balances in (10 a,b) as well as contact line 
condition (|13j) arise from conservation of momentum. The remaining surface equations are 
all derived from requiring nonnegative entropy production in the surface phases, with results 
from mass and momentum conservation also applied in certain circumstances. Whilst this 
connection to entropy production is not immediately apparent from the forms of the equations, 
the derivation connects terms involving the entropy to the surface chemical potential p s , 
arising from the Gibbs relation. The interface formation model equations (7b) and (10c), and 
consequently (8a) and (11a) are then obtained through the assumption that the deviations of 
the surface density from equilibrium are small so that the chemical potential may be expanded 
in a Taylor series about p s e , the equilibrium surface density (see §4.3.2 of Ref . y, or equation 
(21) of Ref. |), and thus 

dp s 



= + dps 



(p S -p S e), (15) 

P"=P 3 e 

so that the governing equations are now written in terms of the surface density — with 
dpsp s (pl) being absorbed into a surface coefficient becoming the relaxation timescale r. It is 
this treatment that then motivates a similar consideration of the line phase of Bedeaux by 
Billingham, 23 and why although condition (|14|) may initially appear more like a momentum 
balance, it does in fact arise from energy considerations, effectively suggesting that the flow 
through the contact line is driven by a difference in the chemical potentials of the interfaces 
there. We reiterate that this additional condition requires treating the contact line region 
as a separate thermodynamic entity, which Shikhmurzaev disagrees with, suggesting that 
there are no singularly strong forces when going down a dimension from surface to line in 
an analogous way to the asymmetric intermolecular forces when going from bulk to surface, 
shown in Fig. [IJ which give rise to the surface tension. 6 

Before proceeding with the investigation of the interface formation model equations for the 
droplet spreading problem of interest, we will make some further remarks on the model. There 
are a number of "non-classical" ingredients which are included within the model, having arisen 
through its derivation using non-equilibrium thermodynamics. In particular, these include 
generalised slip (in (7a)), mass transfer between bulk and surface layers (in (7b), (8a), and 
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(10c), (11a)), and flow-induced Marangoni effect (within (7a), (8b), and (10b), (lib)). The 
model, as mentioned previously, will reduce to the Navier-slip model if the relaxation of the 
surface tension is instantaneous. The slip terms remove the moving contact line singularity, 
although the subtle distinction is that (as highlighted in Fig. [2D the slip occurs between bulk 
and surface layers rather than at the solid surface. The mass transfer and Marangoni effect 
terms are then intrinsically linked as it is the finite time surface tension relaxation which 
induces the Marangoni effect, and the surface tension varies dynamically due to the variable 
mass in the surface layers. These features are critical to the ability of the interface formation 
model to allow and determine the dynamic variation in the microscopic contact angle, as it 
evolves with the dynamic surface tensions through the Young equation Q. 

It is certainly of interest to consider whether one surface may be described classically, 
with the other using elements of the interface formation model. The assumption made here 
that the model parameters a, /?, 7 and r are equal for both surface layers would need to 
be relaxed, and further experimental work to determine the relaxation time r for different 
surfaces is required to justify the assumption that it is negligibly small. A previous study by 
Monnier and Witomski— considered the interface formation model numerically for a plunging 
tape. They assumed that the liquid-gas surface tension was always at its equilibrium value, 
so some analysis of a simplified model has been undertaken. The authors there only make this 
assumption for simplicity, however, and state that their results are only a first step (suggesting 
that the model is too simplified to interpret results from a mechanical point of view), with 
development of the full model in progress. 



2 Application to droplet spreading 
2.1 Governing equations 

We now formulate the governing equations for our prototype system, a two-dimensional 
droplet on a horizontal planar substrate. (X, Z)-axes are chosen parallel and perpendicu- 
lar to the planar solid substrate, respectively, with X = the centre line of the droplet, 
X = ±A(T) the location of the contact lines, Z = H(X, T) the height of the droplet, and 
u = (U,V), see Fig. [3) The bulk motion of the incompressible Newtonian fluid is given by 



Z 




Figure 3: Symmetric two-dimensional droplet of thickness Z = H(X, T) with contact lines at 
X = ±A(T) and coordinate axes chosen to be parallel and perpendicular to the planar solid 
substrate. 
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the continuity equation, dxU + dzV = 0, and conservation of momentum equations, 
'dU TT dU xr dU\ dP fd 2 U d 2 U\ 

32 v q2i 



dT 



dX 



fdV dV Tr dV\ dP (d 2 V d 2 V\ 

p {^ + U dx + V dz)=-dz + ' l {dx^ + d^)- 



(16a) 
(16b) 



The solid-liquid interface has unit normal ns = (0, 1) into the liquid, v s s = (u s s ,0) from ([6]), 
and 



dU dV 



+ = 



ldai 



2 dX 



pV = Xr 



dps , _ ^ - P| e 



dT 



dX 



2 u+a di- 



PS - PSe 



v s s = 7(p s i0 )s-Ps)- (17) 



The gas-liquid interface has normal pointing into the liquid droplet 



1 d 2 H 

— (d x H,-l), giving V s -n G = ^|^-, 



where % = ! + (dxH) 2 is written for simplicity, and then 



n G • T • n G = ^ 



dU 



so that the gas surface equations become 

dH 



dX 



dH 



dU dV\ dH 



{%- 



dZ dX J dX 



&f + UG dx 



(n 



P G - P + n G ■ T • n G = a s G V s • n G , 

'— + —\ +. A—— )ffj<^ 



+ 4- 



G 

H dX ' 



(U-u s G ) 



dZ dX J dXdX 

dH v + v G = Xm Vn p -c^ 

PT 



dX 



dT U 



4/3 



1 + 4a/3 



d , „ . . dH d , „ „ . 
uh-U+{v s G -V) — 



PGe ~ PG 



da, 



G 



dX ' 



a, 



G 



7(P 



(0)G 



(18) 

(19) 

(20a) 
(20b) 

(20c) 

(20d) 

(20e) 

(20f) 
(20g) 



where v G = (u G ,v G ). The surface gradient and surface divergence from © have been used, 
and we confirm n G • V s a = 0, n G • Vfn G = 0, and n G • V,v G = 0, as would be expected. 
Finally at the contact line 



Pg [«r - U cl) cos 9 d + v s G sin d ] + p s s {u% - U C l) = 0, 
a s G cos 9 d + a% = a s Ge cos 6» s + a s Se , 
PH u s ~ u cl) = U (p s s - p s Se + p s G - p s Ge ), 



(21) 
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as Ucl = (Ucl,®)- Noting that the free surface moves at velocity -v^ • ng, then we also 
have: 

• n<7 v s n — u s n dxH , . 

Ucl = = . Q G it ■ ^ 

sin 6 d sin 6 d V7-L 

It is possible to simplify the above equations by removing excess variables, details are recorded 
in the Appendix. The next step will be to nondimensionalise these equations and determine 
which terms can be neglected under the long-wave approximation. We note that for our 
specific droplet spreading problem we have conservation of mass suggesting 



q f A(T) f H(X,T) f A(T) 

Of J o j o pdZdX = ^j o (Ps ~ Pk + PG- PGe)dX, (23) 



taking into account the mass flux into the surface layers. There will also be a number of 
additional constraints imposed to preserve symmetry, relate microscopic contact angle to the 
slope of the free surface etc., but these will be described later. 

2.2 Nondimensionalization 

We proceed by defining nondimensional variables using the following substitutions: 

X = Lx, Z = eLz, T=~t, H = eLh, U = Uu, V = eUv, U C l=UU C l, 

LA 



P=^P, Ps = PSe + ^Ps, Pb = Pbe + ^Pb (2.1,, 



where lower-case letters and barred variables are nondimensional quantities. L is a typical 
length scale in the x-direction such that x = ±A/L = ±a is the location of the contact lines, 
U is a typical horizontal velocity scale, and e is a small parameter allowing for the long-wave 



approximation (often referred to as the "film parameter" in thin-film studies, e.g. Ref. [60j). 
Here we assume e = 9 S , with Od = 0(e). Our governing equations contain the following 
dimensionless numbers 

pUL 7 x 10 7 Le 3 - /3Le 10 10 Le - 3UL 3 x 1(T 4 

Rg — ~ , b — ~ , /in — — ~ , 

fx 3 ' H 3p 3 ' 7rRee 2 2e 2 

1 _ 3tU\ 2 x 10~ 6 e 

a = ap 



12' Le 3 3L ' 

> = > = £^ °" = t — °-- <*> 

where certain scalings have been chosen to achieve the fullest balance in the equation for the 
pressure jump at the gas surface, and U = a s Ge e^ / (3p) = O(70e 3 /3). The nondimensional- 
ization is motivated by that of Billingham 2 — with additional scalings based on Hocking— to 
allow for direct comparison to the Navier-slip case, and hence why certain numerical values 
in (|25|) are present. Approximate parameter values are based on those given in Table 1 of 
Billingham— for water at room temperature, and relevant quantities are also included in 
Table Q] of the present communication for self-consistency. 
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Table 1: Physical parameters estimated for water at room temperature, values from 
Billingham, Ref. 1 24. 



Physical quantity 


Symbol 


Estimated typical value 


inverse surface layer compressibility 


7 


2 x 10° m 2 s~ 2 


surface layer relaxation timescale 


r 


1(T 8 s 


bulk fluid density 


P 


1000 kg m" 3 


equilibrium surface layer density 


PSei PGe 


10~ 7 kg uT 2 


bulk fluid viscosity 


P 


10~ 3 kg m^s -1 


coefficient of sliding friction 


P 


10 7 kg urV 1 


Darcy- type/sliding friction coefficients 


a/3 


1/12 (dimensionless) a 


entropy production constants 


Uq, U G0 


10~ 4 m s" 1 


equilibrium surface tension 




7 x 10~ 2 N m" 1 



a This value is obtained by Shikhmurzaev 4 " through an analogy with flow in a 
plane channel for the shear flow in the surface layer, and is in turn used by 
Billingham.— >2Z 



Considering e <C 1, we have at leading order in the bulk 



du dv dp d 2 u dp 

dx dz dx dz 2 ' dz 



provided f Re<l, and at z = on the solid surface 



du 1 dp 



5 



dz 2 dx 



3/3 u, 



XmRo 



-Ps, 



a 



d 



2 t:S 



3 d 



f3 dx 2 2 dx 



Ps 



(u\z=o) 



Ph 



(26) 



(27) 



making the assumption e 3 <C Xp. Then on the gas surface, at z = h, at leading order 



dh dh XmRo s 

m +u dx- =v + — pG 

32 



-P 



T 



6Ph 

dx 2 ' 

l + Aad 2 p s G d_ 

4^ dx 2 6 dx Wz=h \ 



du dp s G 
dz dx ' 



PG, 



(28) 



assuming e 3 <C A, e 2 Ro <C 1. Finally at x = a, z = 0, the contact line conditions are 



1 + 4a dp s G - v / u a dp s s 



U CL 



0, 



(Pb + Ps) = l 



dh\- 

dx J 



P 



u 



a dp° s 



Ucl 



Cao 



(Ps + Pg) 



2 3/3 dx J 3A 

where the same assumptions as at the solid and gas surfaces have been made, and 



Ucl = u 



dh 

dx 



. XmRo -s 



(29) 



(30) 



The typical values given in (|25p support our assumptions of e 3 <C Xp, e 3 <C A, e 2 Ro -C 1, and 
e 2 Re <C 1. The inclusion of mass transfer throughout the equations of the modern interface 
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formation plays a role in the contact line velocity, as seen above. This can be thought of 
physically as being caused by the interface between the gas surface and the bulk moving at a 
different velocity to the actual droplet interface due to this mass transfer. In particular, we 
have conservation of mass from (|23|) suggesting 

| jT hdx = J\p% + ph)dx. (31) 

For comparison, we record the form of the governing equations if Navier-slip with 9^ = 9 S is 
assumed. In this case the bulk equations are as in ([26]) . and at the surfaces we would instead 
require 

/3ns du , du d 2 h dh dh 

at z = : v = 0, u=—- — , at z = h : — = 0, -p = 3—— — + u — = v 

3 oz oz ox A at ox 

(32) 

which may appear if the surface layer densities are equal to their equilibrium values, i.e. 
p% = Pq = 0, and (3ns = P ■ The reduction of the conditions at the contact line for Navier- 
slip is less obvious, but if in addition p = (so that either the fluid density at the solid surface 
is zero, or infinite at the gas surface) they reduce to u = Ucl and 6^ = S . 

It is of interest to note the regime when the additional terms in the Xm = 1 formulation are 
insignificant, occurring when .Ro <C 1. Based on the values in (|25p this is when 1.2 x 10 -2 <C 
e <C 1, so the difference between the two interface formation models is negligible if we take e 
in the above range. 



3 Asymptotic solution of the long-wave equations for Xm = 

We first consider the regime where Xm = separately from Xm = 1, as it is in this original 
formulation of the model where the additional contact line condition has been required in 
the analysis of Billingham.— >2Z The analysis for Xm = 1 follows similarly, with details given 
in Sec. HI The bulk equations ([26]) and the second condition in ([28]) are solved to find 
p = p{t, x) = —3d^h, and 

3z 2 d 3 h nA ,^ . _ , z 3 d*h 3z 2 dA dB 

U = - — Q^ + 3A ^ x > + ( ' ^ V = ^dx^-—^- Z ^ (33) 

where we have removed an additional term arising from an integration through the use of 
v = on the solid surface (z = 0). The quantities 3A(t,x) and B(t,x) are the shear stress 
and the slip velocity on the solid surface respectively. The remaining conditions on the solid 
surface combine to give 

This equation in the Navier-slip model gives B = /3nsA, so could be considered as the limit 
f — > 0. Turning attention to the gas surface equations, on z = h, the first two equations 
in (|28|) for Xm = are the same as those in the Navier-slip case leading to the free-surface 
equation: 

+ + (35) 



dt dx V 2 dx 3 2 



14 



The other density conditions simplify to 



h ^h-A = f— 

dx 3 dx 2 



3 2 d 3 h 
2 8x^ 



3Ah - B + 



l + Aa f, d 3 h 



4/3 



h 



dx 3 



A 



(36) 



having eliminated the p s G dependence. The equivalent Navier-slip condition is A = hd 3 h, 
again found in the limit f — > from the interface formation boundary condition. The above 
considerations for Navier-slip lead to the single equation for the motion of the free surface of 



dh + d_ 
dt dx 



h 2 (h + p NS ) 



d 3 h 
dx 3 



0, 



(37) 



as expected. Finally the requirements at the contact line (x = a, z = 0) are 



1+'-)B-U C l(1 + p) 



dx J 



ZD, 
(38) 



where Ucl = B(a), and here D = (p s s + p s G )/3 satisfies 

m d \\ + Aafd 3 h A -\ 2pa A 3 l2 d 3 h 



3Ah- B 



(39) 



We note here that A is removed in the first contact line condition when 1 + 4a — 8pa = 0, 
a value which is found to suggest an unphysical infinite shear stress when considering inner 
regions near the contact lines. It is expected that either the combination of values causing 
l + Aa — Spa = is also unphysical, or that the model assumptions break down in this case — 
such as the assumption that a, j3 and 7 take the same values at both solid and gas interfaces. 
We continue by taking parameter values away from this specific combination. Considering 
the contact line conditions for the Navier-slip case we take f — > and notice that the second 
equation yields (d x h) 2 ~ 1, suggesting that the microscopic dynamic contact angle is equal 



to the static contact angle (6 S = 6^) at the droplet rim x 
line conditions are automatically satisfied if p = 0. 
We then also require the conditions 



±a. The remaining two contact 



h{a) = 0, 



dh 

dx 



0, A(0) = B(0) = 0, 



x=0 



f 

Jo 



hdx 



(40) 



corresponding respectively to zero free surface height at the droplet rim, a symmetry condition 
for the free surface at the droplet centre, zero shear stress and slip velocity at the centre of 
the wall (again by symmetry arguments), and the requirement that the cross-sectional area 
of the droplet remain constant. 



3.1 The outer region, (3 3> 1 

For convenience, we transform the domain —a < x < a to —1 < y < 1 via the mapping x = ay, 
and as such we are now able to write the governing equations in terms of the (non-moving) 
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coordinate y, as 



dh 


ay dh 


i d 


' h 3 d 3 h 


dt 


a dy 


a dy 


2a 3 dy 3 



j3B 

h d 3 h 
a 3 dy 3 



A 
A 



rp d 2 
a 2 dy 2 

a 2 dy 2 



1 + 4a 
4 

3h 2 d 3 h 
2a 3 dy 3 



- -Ah 2 - 
2 

b-Ia 

(3 . 



Bh 



3Ah -B + 



1 + 4a / h d 3 h 



4/3 \a 3 dy 3 



A 



(41a) 
(41b) 
(41c) 



where a 
line (y = 

2a 



da/dt corresponds to the spreading rate of the contact line, with at the contact 



i). 



jWB-A)-- 
where 

_ f d 

a dy 



B 1 + 4a 



4/5/3 



h d 3 h 
a 3 dy 3 



A 



2V = 1- 



1 dh 

a dy 



1 2 



(42) 



1 + 4a ( hd 3 h 



4/3 V« <9y 3 



A - 2,0/35 



2pa , 3d 2 <9 3 /i 
+ -^^4 + ^7^7 - 3Ah 



2a 3 dy 3 



B 



(43) 



Slip is not significant in the outer region, and as such we take the limit /3 — > oo (with all 
calculations then at leading order in an asymptotic expansion in /3 _1 ), and at the same time 
consider the quasistatic limit |d| <C 1. We introduce quasistatic expansions of the form 



a(t) = a (t) + ai(t) H , 

h(x,t) = h (y,a ) + hi(y,a ,a ) + h 2 {y,a ,d ,ai,dl,a ) H 

A(x,t) = A (y,a ) + Ai(y,a ,a ) + A 2 (y, a , a , Oi, d§, a ) + • ■ 
B(x,t) = j3~ l [B Q (y,a ) + 5i(y,a ,d ) + B 2 (y, a , d , ai, dg, 



ao) + 



(44a) 
(44b) 
(44c) 
(44d) 



where the time dependence of /i, A, and 1?, enters through a and its time derivatives. As 
written, we have assumed that corrections to the leading order radius satisfy |ai| <C |do| "C oo> 
implicitly assumed by Hocking, 51 and Savva and Kalliadasis. 61 

For clarity here, we state the assumptions that will be required in our asymptotic analysis 
and outline the procedure to follow. There will be two main asymptotic regions. An outer 
region in the bulk of the droplet, where slip is not significant and an inner region at 0(1/ (3) 
near the contact line (only one contact line being investigated due to the symmetry of the 
droplet). An intermediate region matching these two was found necessary for the Navier- 
slip model by Hocking,— but here we are able to match the variables through their cubes, 
or negative cubes (depending on the variable at hand), with an intermediate region merely 
justifying this procedure. A quasistatic expansion is then performed in both inner and outer 
regions, and to extract the information about the spreading rate, we will find that only ao 
and do will appear. The assumptions explicitly made in the procedure are f = O^^ 1 ), 
| a-i | <C |do| -C ao, /3 -1 "C |do|, do = 0(oq), and will be commented on further where they 
arise. The final assumption listed is merely due to the way in which the quasistatic expansions 
are applied, with terms involving a±, do and d§ occurring at at higher order than 0(do), so 
that a\ = 0(do) = 0(oq). As we are only interested in the leading order terms in j3, all of the 
corrections to the radius and velocities are all greater than 0(1/(5) as written. This is found 
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to be consistent, as do depends only logarithmically on j3 (as originally found by Hocking— 
for Navier-slip). 

From (|41bp we see that B = o(l), motivating the scaling of B with as given in 
the expansion above to achieve the fullest balance. We make the further assumption that 
f < 1, which is reasonable based on the parameter values (|25|) . This will be confirmed 
after considering the inner region, where we find it necessary to have r = 0(/3 _1 ) to achieve 
d x h\ y= \ = 0(1). These considerations yield at leading order 



Ba 



hp s hp 
an dy 3 



dy ^ "° dy 3 



l\h 



0. 



(45) 



with conditions 



A (1) = 0, 



Oho 
dy 



ho(l) 



0, ^ 



y=l 



dho 
dy 



y=0 







1 1 

hpdy = —. 
ap 



(46) 

Equations ()45p and conditions (|46p (apart from the contact angle condition) are solved to 
find 



ho 



3 1-1/' 



2a 

The equations for the next term are 
ho d 3 h\ 



giving 



An = Bn 



0. 



(47) 



A 1 =B 1 

d 3 h x 
dy 3 9(1 — y 



0% dy 3 
4d ajjy 



with solutions 



2\2 



A x = 

dhi 
dy 



B 



2a a y 



1 ■ 5 

g a oa 



3(1 -y 2 ) 



In 



1 + y 



1 



y 



3y 



(48) 



having used dthn = and ao ho, and where arbitrary constants have been determined from 
h\(l) = and (d y hi)(0) = 0. The solutions for A% and B\ automatically satisfy the symmetry 
conditions -Ai(O) = ^i(O) = 0. The behaviours of the free surface slope, A\ and B\ as we 
approach the contact line at this order are then found as 



dh 3 1 4l 



e 3 (a 



2a 



A ~ -a al(a 



x 



B 



3/3 



aoal(a - x) 



(49) 

as x — > a. The behaviours of A and B in (|49p include only one term in the expansion for 
|d| < 1 as the Oil) terms are zero. As is seen in a number of other works with slip 51,61 and 
precursor films*" for the matching of d y h, the logarithmic terms appearing in the second term 
must be included for the correct matching procedure. We expect the same to be necessary 
for the matching of A and B, so we continue to O(|do| 2 ) to find terms of that size. Returning 
to (|4ip with expansions (|44p . we have 







4 



Ao=Bo = ^r 



dh\ dh\ 

a o^ 1" a o^^ 

aao dao 

, d 3 h 2 d 3 h, 

ho^r^r + h l 



o dhi d 
a ° aoy ^y- + dy- 



hn 



..d 3 h 2 
dy 3 



+ 3h^h 



d 3 h x 
dy 3 



dy 3 



dy 3 



(50a) 
(50b) 



having used (|48p . We may have expected terms involving the correction to the radius a\ in 
the above expressions, but it is readily seen that the terms arising are O(aido), and hence do 
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not feature at this order. Using ([IT]). (jIHj) and applying J* h\dy = 0, these may be solved to 
determine the behaviours 



7a2 



Ao = Bo 



aka. 



r e 2 (l-y) 



3 4 (l-y) 



In 



as y 



(51) 



noting that there are no terms involving ao at this order. In the expansions for A2 and B2, 
they occur first at O(do ln(l — y)), and are smaller than the terms kept due to the assumption 
that Sq = O(clq). Collecting all the results together, our two term matching conditions become 



1 and 



8a2 



3 a 



3 4 a 



In 



x 



dh 3 1 . 

7^~72 +n a 0«0 ln 



Ox 
e 2 (a 



'0 

x) 



9 



e 3 (a 



2a 



B 



2a 
d ao 



1 and. 



8a2 



0"0 



3/3 ao — x 3 4 /3 ao 



In 



e 2 (a 



2a 



(52) 



as x —7- a. These asymptotic expansions will be matched to their counterparts from an inner 
region close to the contact line to determine the droplet spreading rate. The matching is 
carried out in a region such that for A and B the O(do) terms occur at higher order than the 
0(clq m(ao — x)(clq — x) _1 ) terms and the O(do(ao — a;) -1 ) terms respectively, and equivalent 
restrictions for d x h, to ensure asymptotic validity. 

3.2 0(1/(5) Inner regions 

The outer solution found does not satisfy the boundary conditions at the contact line, and 
as such motivates the presence of a boundary layer, or inner region, close to the rim of the 
droplet. In fact, for the slope of the free surface defining the contact angle to satisfy the 
boundary conditions, we would require a 2 = 3. If the initial droplet radius satisfies this value, 
the droplet will remain in equilibrium, but radii different from this will cause motion of the 
free surface to return to this equilibrium. Considering distances a — x = 0(1/ P), we make 
the change of variables 



x = a — £/3 



h = f/r 



(3A, B = B, 



(53) 



where the scalings for h, A and B are motivated by their forms in the outer solution. We 
then find the governing equations become 



1 a* . _ d 



1 o9 3 $ 3 



de+ -A^ + m 



~ ~ - d 2 

A - B = (3f Pde 

J 3 $ , 5 d 2 
^ h A = Br 



aA 



l + 4a 



B 



3 t2 <9 3 1' * t a l + 4<5 



d 3 ^ j 



(54a) 
(54b) 
(54c) 



with, at £ = 0, the contact line conditions 
l + 4a 



B + 4a(A-B) 



2p 



a 3 ^ j' 



2V = l-[^) , 4a(B-A)-B 



2 Cap 



V. 



(55) 
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where 



V = f3f 



d_ 



+ >4 + 2pB - 2paA + — + 3Atf + 5 



2 



(56) 



We then also require ^(0) = and that the solution matches to the outer region. The 
dependence on the combined parameter (3f in the above equations is clear, and that to 
construct an asymptotic solution with the contact angle 0(1) requires f3f = 0(1) as assumed 
in the outer region. Noting that fif = 0(2 X 10 3 e 2 ) from ([25]) for water at room temperature, 
the above assumption is a reasonable one, but we will consider a range of values. It is 
noteworthy that /3f was also found to be an important parameter in the gravity driven thin- 
film flow situation.™ As seen for the Navier-slip model,— and also in the comparison between 
slip and precursor film models,— the inner region solution as £ — > oo will determine the 
matching condition to the outer region and provide the relationship between the spreading 
rate d and the radius of the droplet a. We will proceed by considering a quasistatic expansion 
for the inner equations, but first we note that there is a correspondence between the limit 
£ — > oo and the situation when j3f «C 1. In particular, given that A and B cannot display 
exponential growth as £ — > oo, and also that for matching to the outer region solutions (|52p 
we require Vl//£ 2 — > 0, A — > 0, and B — > as £ — > oo, then we have to leading order 



A ~ B 



d 3 ^ 



i a* .a* 



d_ 



tf 2 (l + * 



a 3 ^ 



as £ 



oo. 



(57) 



3.2.1 The inner region with quasistatic spreading 

To progress with the analysis of the inner region governing equations ([31]) . we assert that 
A = O(do), -B = O(do). This is motivated by observing that ^4 = O(do), -B = O(do) in the 
outer solution, thus for solutions of A and -B in the inner region to match, they must be at 
most O(do). Alternatively, the original variables A and B correspond to the shear stress and 
velocity at the solid surface so would not be expected to contribute to the leading order static 
terms. Introducing the quasistatic expansions 



a(t) = a (t) + a\(t) H , 

t) = #o(£, a ) + d , oo) + *2(£> OO) oo, oi, do, do) H 

i(x,t) = ii(£,d ,a ) +A 2 (€,ao,ao,.. •) H , 

B(a;,i) = Si(£,d ,a ) + -B 2 (£, a , o , ■ . ■) H , 



(58a) 
(58b) 
(58c) 
(58d) 



then 



- ([3U]) at leading order are 







_9 fI^ 3 5 3 *o 



9£ V 2 ^ U ^ 3 
with the contact line conditions 

9 3 * 



<9 3 ^c 



ae 2 



o, 



2 ^4 



9 3 ^o 



a 3 ^ c 



3 T 1 + 4a 
2* + — 



9 3 ^o 
' d£ 3 



3 1 + 4a 
2* + ^ 



3 1 + 4a 
-^0 + : 



d^ 



(59) 



(60) 
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at £ = 0. As ^o(O) = 0, we find the first equation of fl59|) gives # = K £ 2 + K^, which 
automatically satisfies the second equation and the first and third contact line conditions. 
Requiring ^o/^ 2 - * as £ — > oo to match to the outer region and applying the second contact 
line condition gives = £• We note that the third contact line condition in (|6U|) is that 
added by Bedeaux— and Billingham,— >2Z and whilst it is automatically satisfied here, it was 
not required to determine the solution at this order. 

Having determined the leading order free surface behaviour in the inner region, we consider 
the next order in the governing equations, where they satisfy 



At 



a 



d_ 



l^d 3 ^! 3 j 



2^ ae 



a.A\ 



l + 4a 



-A 



Bi 



+ 3Ai^ + B 1 + 



1 + 4a 



^ 3 ^i A 



(61a) 
(61b) 
(61c) 



having neglected the f3 1 dt i &o = do£(/3ao) 1 term as |do| <C /3|do|. At £ = 0, we have 



B 1 +Aa{A 1 -B 1 



1 + 4a 
2p 



^ 3 *i i 



4q(5i-Ai)-Si 



2 Cap 



V. 



where 



V = /3f 



l + 4a , - 



2paAi + -|— — ^ + 3Ai^ + Bt 



(62) 



(63) 



with *i(0) = 0. We also require ^i/^ 2 — > 0, Ai ->■ 0, and i?i — ► as £ — ► oo to match. As 
noted previously, it is the behaviour as £ — )• oo that is crucial to the leading order spreading 
dynamics of the droplet, so we proceed by considering an asymptotic expansion of the form 



d 



i (c i iln£ + Ci 2 ), 



8=0 



with similar expressions for A\, B\, to find 
*i = d [£ ln£ + c £ + 0(ln£)] , A x = d [r 1 + 0(r 2 )] 



(64) 



= «o [r 1 + o(r 2 )] , 

(65) 



as £ — >• oo, in agreement with (|57|) at O(do). The constant Co is unknown here — it will depend 
on all of the parameters of the problem, and must be determined from the full solution of 
(|6ip - (|63|) . considered in Sec. 13.2.21 The required two terms have been found for the matching 
of the slope in the inner region. As in the outer region, we require terms in both do and clq 
(and possibly other terms at this order) for A and B to also have a two term expansion. From 
(foT|). we have that 



An ~ Bo 



d 3 ^ 2 



of 



d_ 



d£ 3 



(66) 



as £ — > oo, where we have neglected the terms including by making the assumption j3 1 <C 
| do |, an assumption also made by Hocking.™ This assumption means that our asymptotic 
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analysis is not valid very close to equilibrium, although we will see from our comparison with 
numerical results for the full problem that excellent agreement is still found. Using (|65|) and 
combining all of these results in terms of outer variables 

-|^~l + a ln[e 1+c ^(a -x)] , A ~ $B ~ ln[e c ^(a - x% (67) 

ox L J ao — x ao — x 

as /3(a — x) — > oo, which must be matched to ([52]) . It is apparent, as for the Navier-slip 
situation, that the logarithmic terms of the free surface slope (and of A and B) in the inner and 
outer regions do not generally match. This may be resolved by considering an intermediate 
region as was done first by Hocking, 51 however Savva and Kalliadasis noted that matching 
may be resolved by considering (d x h) 3 between inner and outer regions, with the intermediate 
region simply justifying this step.— Equivalently, here an intermediate region motivated by 
the analysis of Hocking merely justifies the matching of (d x h) 3 , A~ 3 and B~ s . 

A comparison of the cube of the free surface slope and the negative cube of A and B from 
the outer expansions (|52|) and the inner expansions (|67p then gives the relationship between 
the spreading rate do and the droplet radius ao as 

/ e 2-c \ 

3d a[jln^_J = a i -27. (68) 

This spreading rate has a similar form to the result for Navier-slip, where 0ns = 1/13 and 
Co = 0. This means that the interface formation model (with Xm = 0) has a spreading rate 
equivalent to a Navier-slip model (with 9d=0 s ) with slip length (3ns = e _c ° /(3, although we will 
also see an equivalence with other slip models with velocity dependent microscopic dynamic 
contact angles. From (f68j) it is evident that provided ao is away from equilibrium, a(oo) = \/3, 
then we have ao = (ci) 1//7 , where c = — 63/ ln[e 2_c °/(2ao/3)]- Neglecting the weak logarithmic 
dependence of time on c, this gives agreement with the spreading of Tanner's law. 62 1 63 As the 
equilibrium radius is approached the behaviour becomes exponential, so that ao ~ a 00 — le~ wt , 
where e <C 1 is some constant and 



2 

w = — 

V3 



2vW 



(69) 



These spreading behaviours of the droplet radius both reduce to equivalent expressions for 
Navier-slip when Co = 0, as expected (see Savva and Kalliadasis, 61 taking the planar substrate 
situation). Figure IHa) shows the evolution of the droplet radius a(t) for various values of cq 
from initial radius a(0) = 1 and j3 = 10 3 using the spreading rate found in (|68p . showing that 
increasing cq delays the spreading process, but that the final equilibrium radius is unaffected. 
The asymptotic expansions (|67p are valid for 

x » a — i exp ( — !— — 1 — cq ) , (70) 



P Vl a o 

and for matching between inner and outer regions we have £ = f3(a — x) S> 1. Together these 
require 

(l + c )|d | « 1. (71) 

This reduces to the familiar requirement |do| <C 1 for Navier-slip, Co = 0, but imposes a greater 
restriction on the spreading rate for positive cq. Using (j68p . we show the critical value of the 
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initial radius a c (0) in Fig. Hl^b), so that when droplet spreading starts closer to equilibrium 
than that value, the asymptotic procedure is valid. As noted above, we also required /3 _1 <C 
| do | for the asymptotic analysis to hold. Once the droplet has spread sufficiently for this to 
be invalidated, the droplet will already be negligibly close to equilibrium, as our analysis is 
based on the assumption that ^ > 1. We note that an asymptotic investigation to relax 
the quasistatic assumption has been performed by Eggers for a slip model in the vicinity 
of a contact line.~ Even for this simpler model, matching to an outer flow region (such as 
required for droplet spreading) was deemed a challenging and open problem, and thus it is 
also beyond the scope of the present work with the interface formation model. 
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Figure 4: (a) Evolution of the droplet radius a(t) for various values of cq from initial radius 
a(0) = 1 and j3 = 10 3 showing delayed spreading for increasing cq. (b) The critical initial 
value of the droplet radius a c (0) for varying Co and f3. Asymptotic validity if a(0) is closer to 
the equilibrium radius, a(oo) = \/3, than the value of a c (0) here. 



3.2.2 The full O(do) inner region behaviour 



We drop the subscript from ao and do for the remainder of this section as we have found 
that corrections to the radius do not enter into our analysis. Returning to the O(d) inner 
region equations (fol~j) ~ ([63]) . we now wish to find a full solution to complete our asymptotic 
description and obtain an 0(d) correction to the microscopic contact angle. Equation ()61ap 
may be solved to give 



Si 



1 2 dH l 3 



-X 



d£ 3 



AX, 



(72) 



having used — > 0, — > 0, BX — > as £ — > 0, to remove a constant of integration — the 
latter two conditions being imposed as the shear stress and velocity (or correspondingly A 
and B) at the solid surface should at worst be logarithmically singular as the contact line is 
approached (they will, in fact, be found to be finite in later calculations). This determines 
-Bi(O) = d, and reduces the problem to a seventh order system of two ODEs 



(2 + 36ii + e 2 ^ 



2d + 



Pfp d 2 



4 



[I + 4a)e 



+ 



3C 



1 + 4a 



+ 8a) A 



, Pf d 2 



(4f + 1 + 4a) + (6£ + 1 + 4a) A x 



(73a) 
(73b) 
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subject to at £ = 



1 + 4a - Spa t . d 3 *i 2p(l - 4a) . 



1 + 4a 



-V 

(1 - 4a)a + 4a^i 



d£ 3 1 + 4a 



-a, 



2Ca 



P. 



with *i = 0, and where 



_ /3rd 

2? = 

4 9£ 



(6£ + (1 + 4a) (1 - 3£p) - Spa)A 1 + (4£ + (1 + 4a)(l - £p))£ 



<9 3 *i 



(74a) 
(74b) 
(74c) 

(75) 



We note that the boundary conditions are able to support a solution with (£<9| , I , i)(0) = 0, 
which allows the model to alleviate the logarithmic pressure singularity associated with the 
Navier-slip model. This weak singularity occurs at the contact line for Navier-slip (seen 
for example, in Refs. [7 65), where a finite force is nevertheless predicted. We will return 
to this issue when considering the asymptotic behaviour of the interface formation model 
as we approach the contact line, where comparisons between models will be clearer. The 
ODE system (|73p is subject to the contact line conditions at the singular point £ = and 
the matching conditions as £ — > oo of A\ — > 0, B\ — > 0, *i/£ 2 — > 0, all being satisfied by 
behaviour ([65]). We perform a local analysis about these points separately to determine the 
relevant asymptotic behaviours. 



Eigenmode analysis: We are interested in imposing the contact line and matching be- 
haviours 



i ~ Pot A x ~ So, as £ -»• 0, 



and 



d [ln£ + Co] £, A\ ~ d£ , as £ — > oo, 

(76) 

on the ODE system (|73p , where Pq and So are the coefficients of the leading order terms in the 
expansions as £ — > of *i and A\ respectively. To determine the correct specification of the 
inner problem the degrees of freedom contained within the asymptotic behaviours must be 
found, determining the number of boundary conditions imposed. We perform an eigenmode 
analysis for the behaviour as £ — > oo by considering 



i = a [ln£ + c ] £ + 5^(0, A 1 = dC 1 + ^i(£) 



(77) 



where 5 is a small artificial gauge. Keeping terms of 0(5), the resulting equations in ^\,Ai 
have the seven asymptotic behaviours for the eigenmodes of 
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• (78) 



The second mode corresponds to changes in the parameter co, which has to be determined 
numerically, with the first mode corresponding to another degree of freedom at higher order 
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in the asymptotic expansion. The third mode and the two positive exponential modes are 
inconsistent with the behaviour as £ — > oo, and as such impose three conditions on (|73p . The 
two negative exponential modes are exponentially small as £ — > oo and correspond to two 
further degrees of freedom. This implies that the contact line expansions are not analytic, 
with the presence of exponentially small terms. These expansions take the forms 



H + c 2 exp 
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(79a) 
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(79b) 



as £ — >• oo, where we see the four degrees of freedom, with the asymptotic behaviour imposing 
three conditions on (|73p . For the contact line behaviour, in a similar manner we consider 



vI/ 1 = Po£ + <5\I/ 1 (£), A l = So + 6A 1 (Z), 
and determine the seven eigenmodes as 
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(81) 



The first mode is inconsistent with the asymptotic behaviour and imposes one condition, the 
other six modes are all consistent and all have free parameters associated with them. With 
the three conditions from the large-£ asymptotic behaviour, we now have four conditions, and 
require three more. Three of the six free parameters at the contact line must thus be fixed by 
the remaining boundary conditions. The correct specification of the problem requires further 
terms in the expansion as £ — > 0. Proceeding systematically, we find 



*i ~ p £ + Ae 2 in £ + P 2 e + Ac 3 + A£ 4 , 

's -a 3(l + 4a)(2 J P 3 + 5i) 



So + 5i£ + 



2pf3ra 



8a 



(82a) 
(82b) 



where 
(1 + 4a)P 4 



Pi 16a 2 -24a + l 5 (l + 4a)d 1 + 4a(l - g) ~ (1 - 4a) 2 > 



6/3f 32a 48/5/3fa 48p/3fa 64a 

is used for simplicity and we see the six free parameters (Pq, Pi, P2, P3, Sq, Si). 
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Contact line boundary conditions: We are now in a position to discuss the boundary 
conditions at the contact line for mathematical consistency. Firstly, we consider the two 
agreed contact line conditions (|74ap and (|74bj) . arising from the mass balance through the 
contact line and the Young equation. Condition (|74a|) gives the shear stress at the contact 
line 

2 |d(l-4a)p- (l + 4a)A 



Sn 



1 + 4a — 8pa 
and condition (j74bj) determines 



(84) 



6(1 + 4a)P 3 = 3[p(l + 4a) - 2]S - (1 + 4a - 8pa)Si - 4=h + 2[p(l + 4a) - 4]P X . (85) 

There remains one further free parameter to fix for a correctly determined system. Considering 
the pressure, which corresponds to the second derivative of the droplet height (as seen earlier 
in ([28]) ). we have 

^ ~ 2Pi ln£ + 3P a + 2P 2 , as f -> 0, (86) 

so that if finite pressure is imposed at the contact line, then this would force Pi = 0, which 
fixes a free parameter. This behaviour then determines (^5|^i)(0) = as suggested earlier. 

As an alternative if finite pressure is not imposed, then the additional condition of Be- 
deaux 23 and Billingham ™* (|74cp requires 

- [4aS + (1 - Aa)h]p\ 

P ° = 2^ ' (87) 

fixing Pq instead of P\. Provided one of these two options is chosen, we will have a well-posed 
problem, whereas fixing both would then suggest that we would have an over-determined sys- 
tem. It is important to make clear that the logarithmic pressure singularity associated with 
Pi 0, both for the interface formation model and in an equivalent contact line expansion 
for the Navier-slip model (discussed later in this section) is not a fatal flaw as it is an inte- 
grable singularity, and thus leads to a finite force at the contact line (rather than the infinite 
force associated with applying the classical no-slip condition). Whether the thermodynamic 
consideration of the contact line leading to the additional condition may be discounted in 
preference of finite pressure remains a debate, and was discussed in Sec. [TJ 

We conclude that the equations of the interface formation model are well-posed for the 
problem considered here, by taking one of the two above options. We will briefly analyse 
the alternative of using the additional condition, but consider our droplet spreading problem 
predominantly with the requirement of non-singular pressure at the contact line, which fur- 
nishes the problem with a (physical) condition to make it fully determined. This is chosen as 
a singularity in the pressure at the contact would arise due to the model, rather than from 
physical origins, and by allowing the (admittedly weak, logarithmic) pressure singularity to 
exist would suggest that the model is invalid in the immediate vicinity of the contact line. 
This does however rely on the viewpoint that the additional condition may be dropped or 
discounted (see the discussion in Sec. [T]). 
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Numerical implementation: We consider (|T3|) for the variables Vl/i/d, Ai/a, B\/a so that 
a is removed from the problem, and with fixed values of the parameters a, f3f , and p (and also 
with Cao / A if using the additional condition). This is solved as a two-point boundary value 
problem (BVP) over the truncated interval [£o> Coo], with £o and ^ small and large enough for 
convergence respectively. Imposing the contact line expansions ([82]) supplies four conditions, 
as determined by our eigenmode analysis, with the matching behaviour (|79p supplying the 
three remaining conditions and cq obtained through 

co~ | ^-£ln£j , as£^oc. (88) 

These behaviours will be imposed at the interval end points. The BVP is then solved using 
the MATLAB bvp4c solver which employs a fourth order collocation method based on a three- 
stage Lobatto Ilia implicit Runge-Kutta formula. 66 Figures [SHE] depict typical inner solution 
behaviours for the finite pressure condition, and parameter values a = 1/12, (3f = p = 1. 
Figures [5ja)-(c) show the O(d) corrections to the droplet height, free surface slope, A and B. 
Figure [5jd) shows the behaviour of the solution as £ — > 0, including confirmation that there is 
no logarithmic pressure singularity through showing £9|^i — > as £ — > 0. Figure [6^a) shows 
the behaviour of the solution as £ — > oo, with convergence to the expected large £ behaviour, 
with Fig. 0(b) showing the determination of Co from the and behaviours. The effect 
on Co when varying f3f for selected a and p is given in Fig. [JJ To also demonstrate the effect 
of using the additional condition on the spreading rate through the value of Co, we include 
Fig. [51 which has varying j3f for fixed a and p, and selected CaoA _1 (/3f) 1//2 . We choose to 
vary this combination as Cao an d A appear only as CaoA -1 in the governing equations, and 
combining them in this way with /3f removes the dependence on the long-wave parameter e. 
Based on the nondimensional numbers in (|25|) . we have CaoA- 1 ^) 1 / 2 « 2T 1 / 2 x 1(T 4 . 

Whilst not undertaking an exhaustive analysis of the flow fields in the vicinity of the 
contact line, we do include two situations for comparison, in Fig. [9) Streamlines are shown 
in the frame of reference of the moving contact line, and plotted in outer variables. In (a) we 
show a typical flow situation with parameter values j3 = 10 3 , a = 1/12, p = f3f = 1 where 
the behaviour is as found for slip models, with the fluid flowing down the free surface towards 
the contact line and away along the solid surface. In (b) parameter values are chosen to be 
/3 = 10 3 , a = 1/12, p = 10, f3f = 0.1. This is a situation where the flow-induced "Marangoni 
effect" occurs (see §5.2.3 of the monograph of Shikhmurzaev for details of this effect 6 ), being 
caused by the surface-tension gradient driven by the flow field in the region of the contact 
line. Both flow visualisations use the finite contact line pressure condition. 



A comparison to slip with a velocity dependent contact angle: Matching of the 
interface formation model equations between inner and outer regions has given rise to a 
modified equation for the spreading rate when compared to a Navier-slip model with 64 = 9 s . 
Given that the interface formation model allows the microscopic dynamic contact angle to 
vary with contact line velocity, it is of interest to see how a similar behaviour would change 
the Navier-slip model. As an example, we consider the contact line condition 

(dh \ 2 
— (a) J =l + 2dciV5. (89) 



26 



x 10 




Figure 5: Numerical solution of the boundary layer at O(d) for typical parameter values 
a = 1/12, /3f = p = 1, using the finite contact line pressure condition. Shown are the 0(a) 
corrections to (a) the droplet height, (b) the free surface slope, (c) Ai/a, and (d) the solution 
behaviour as £ — > 0, including confirmation of the absence of logarithmic pressure singularity, 
since £<9f ^i — > as £ — > 0. 



The equations of the Navier-slip model in the inner region give ^?ons = £> an d at 0(a) we 
have 



INS 



[(1 + £) 2 ln(l + £) - £(1 " 2cjv5 + £1*0] > havin g a PP ned 



liVS 



(0) = catso, 
(90) 

and where other arbitrary constants have been found using ^ins/S, 2 as ( -> oo and 
^ijvs(O) = 0. Thus the behaviour as £ — > oo is found to be 
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27, 



(91) 



an equation for the spreading rate. This shows that the variation of 9d with contact line 
velocity impacts the spreading rate in the same way as found for the interface formation 
model, in equation (|68|) . In the Navier-slip case, we can show analytically that the spread 
rate is affected precisely by the velocity dependence of the microscopic contact angle, whereas 
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Figure 6: The behaviour of the solution as £ — > oo for typical parameter values a = 1/12, 
f3f = p = 1. (a) Convergence of the second, third and fourth derivatives of the droplet height 
and of Ai and its first derivative, to their expected large £ behaviour, (b) Determination of 
the constant cq from the and behaviours, both using the finite contact line pressure 

condition. 



the constant cq must be found numerically in the interface formation model. We note that 
the Navier-slip model is not the only slip model that can be reduced to this spreading rate 
equation. The slip model of Ruckenstein and Dunn— yields u = {^ D /h)d z u, on z = in the 
long-wave approximation.—^ This has the inner solution 



T a 



£ln(l + £ 2 ) + (1 - £ 2 ) arctan £ + |(tt£ - 2 + Ac RD ) 



(92) 



where we have imposed d^iRD(0) = crd®, and we find that (f9Tj) holds with cns = crd- 
A comparison between slip and precursor film models is given by Savva and Kalliadasis,— 
detailing a wider class of models which also reduce to equivalent spreading dynamics. We 
note that the Navier-slip model still has logarithmically singular pressure at the contact line, 
whereas the Ruckenstein and Dunn model has finite pressure since 

—^2 aln£ + < + 0(£ ), — 2--af + 0(£)> as £ ^ 0. (93) 

The interface formation model is also able to alleviate the pressure singularity as discussed 
earlier. Another feature of interest is the contact line velocity. For both slip models and for 
the interface formation model with Xm = we have u\ x=a = a. This suggests that when 
considering a frame of reference moving with the contact line, the contact line is a stagnation 
point, and rolling motion does not occur. We next consider the inclusion of the extra terms 
when Xm = 1; and will find that this stagnation point is removed. This then allows for 
the rolling motion seen experimentally as discussed in the introduction, with an illustrative 
schematic in Fig. [10j 



4 Inclusion of additional mass transfer terms, Xm = 1 

Returning to the nondimensional equations under the long- wave assumption (|26p ^ (|3U|) . we 
now consider the case where Xm = 1 for the modern version of the interface formation 
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Figure 7: The effect of varying (3t on cq for (a) varying a with p = 1 fixed, and (b) varying 
p with a = 1/12 fixed, both using the finite contact line pressure condition. 



model of Shikhmurzaev.— The parameter R$ ~ 1-5 x 10~ 4 e~ 2 , based on the sizes of the 
nondimensional numbers given earlier in (|25p for water at room temperature, so that Rq 
( )( I : if e = 0(1.2 x 1CT 2 ). Analysis considering Rq = 0(1) should thus also be performed. 

The differences to the previous analysis for \ m = appear only in the equation for v on 
the solid surface and in the kinematic condition on the gas surface. In a similar manner to 
the Xm = case, we thus determine 

3z 2 d 3 h „ „ z 3 d*h 3z 2 dA dB R , . t . 

" = --8^ + 3A > + B - 1, = T&3--&- 2 lfe + T^ (94) 

and these additions manifest themselves in the free surface equation only, which now becomes 

dh d ( l,,a 3 /i 3 ,,2 \ - m 

where V satisfies (|3^|) . A(t,x) and B(t,x) once again satisfy and ([5U|) and also at the 
contact line the conditions (|38p are again applicable, but with the modified contact line 
velocity Uql = B — Rq {d x h)~ l V. Mass conservation is determined by ([3"T]) . suggesting that 
mass transfer to the surface layers drives the actual contact line velocity from that created 
by the bulk. Given that in the outer region of the droplet f < 1 < ,8, then all additional 
terms are negligible and the outer region solution for Xm = 1 follows as in the Xm = case. 



The inner region: We next consider the inner region with the scalings (|53p and the qua- 
sistatic expansions (|58p . For the leading order solution we find = £j following a similar 
argument as in the Xm = case. Considering the next order in the governing equations, we 
have 

9 A ,3^*1 , 3 , t2 
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Figure 8: The effect of varying (5t on cq when using the additional condition of Bedeaux and 
Billingham, and allowing logarithmically diverging pressure at the contact line, for varying 
CaoA-^f) 1 / 2 with p = 1 and a = 1/12 fixed. 



with (|61b|) - (|61cp . and where T> is given by (|63p . At £ = the contact line conditions are 
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with *i(0) = 0. Once again we require *i/^ 2 — > 0, A\ — > 0, B\ — > as £ — )• oo to match to 
the outer solution. Studying the extra terms in this Xm = 1 situation, it may be seen that 
they will not contribute at leading order as £ — > oo, and an asymptotic expansion confirms 
that the matching results of (|65p still hold, with differences occurring only at higher order. 
Obtaining further terms for the matching of A and B again follows as for Xm = so that the 
matching behaviours are in agreement with (|67p . Given the outer region is in agreement for 
Xm = and Xm = 1, as is the leading order behaviour for £ — > oo from the inner region, then 
the matching and the result for the rate of droplet spreading (|68[) will be the same. 

Returning to the O(do) inner region equations ([9"6"j) - ([9TP with (|61bp - (|61c|) and (|63p. we 
now wish to find a full solution to complete our asymptotic description and obtain an O(do) 
correction to the microscopic contact angle. We drop the subscript from ao and do for the 
remainder of this section as we have found that corrections to the radius do not enter into 
our analysis, as for Xm = 0. Equation ([96]) may be solved to give 
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Figure 9: Streamlines for two flow situations arising from different sets of parameter values, 
in the frame of reference of the moving contact line, (a) shows the expected flow situation 
with no region containing closed streamlines using parameter values f3 = 10 3 , a = 1/12, 
j3f = p = 1. (b) has parameter values chosen as /3 = 10 3 , a = 1/12, p = 10, f3f = 0.1, and 
shows a situation where the flow-induced Marangoni effect occurs. Both flow visualisations 
use the finite contact line pressure condition. 



where K x is a constant of integration which we cannot remove in this case, as we have 

^ = 4^-^(1 + 45-8^)^(0) 
U ' /3fi? (4 + 2p(l + 4a)) ' 1 ' 

which is finite (assuming A\ (0) is finite) as required and as such does not remove the arbitrary 
constant. After substantial simplification, the first and second contact line conditions together 
reduce to 

and m = %^-*2g£$^,<m 

so that Bi(0) = a + R d^i(0). This also confirms U C L = B(0) - R d^i(0) = a, differing 
from the bulk contact line velocity -Bi(O), and allowing rolling motion to take place. We find 
through a similar eigenmode analysis as for Xm = that the third contact line condition (|97cp 
is not required, provided a finite pressure at the contact line is instead imposed. We are left 
with a system of two ODEs of seventh order, which is obtained from (|61b|) - (|61cp . (|63p . and 
(j97|) and substituting in the above results. We will not record this here due to its cumbersome 
size. 

Numerical implementation of this system with the solution for B\ from (|98p follows as 
for the Xm = situation, and we refer to the details of Sec. 13.2.21 Figure fTTI illustrates the 
effect on Co when varying Rofif for selected values of the parameters (/3f,a,p). We choose 
to vary Rofif as based on the nondimensional numbers in f)25|) . we have Roftf ~ 1/3 with 
no dependence on the long-wave parameter e. As Rofif — > 0, cq approaches the value from 
the Xm = case. Figure [12] compares the velocity fields close to the contact line of the 
two formulations of the interface formation model, (a) Xm = and (b) Xm = 1- The same 
parameter values of /3 = 10 3 , a = 1/12, p = 0.1, (3f = 1, are used in both figures with 
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Figure 10: Diagram of rolling motion in the frame of reference of the contact line, with 
interfaces shown schematically as layers of finite thickness. The contact line created by the 
bulk is shown with a circle, where the velocity is shown as non-zero. This occurs for the 
interface formation model with Xm = 1, but not when Xm = or with slip models, where a 
stagnation point occurs. 



additionally -Ro = 1/3 for the Xm = 1 formulation. The predicted low velocity region is seen 
in (a) for Xm = 0, with the contact line a stagnation point. This prevents nanoscale rolling 
motion as a fluid particle on the surface will not reach the contact line in finite time. This 
issue is resolved in the model with Xm = 1, as seen in (b). We note that mass conservation 
for bulk and surface layers is guaranteed from the formulation of the model, as the original 
equations lead to the mass balance (f2"3"|) and the first contact line condition in (|2"Tj) . 



5 Numerical analysis 

Returning to the full problem in both model formulations, we wish to solve the governing 
PDEs numerically to verify our asymptotic analysis. We will first consider the Xm = case, 
with the Xm = 1 case then following a similar procedure and by considering the model with 
the requirement of finite pressure at the contact line. 



5.1 Numerical solution of the full problem for Xm — 

We now return to the governing equations of the full boundary value problem for Xm = in 
terms of the non-moving coordinate y, where x = ay, given in (f4"Tj) - (f4"3"|) . and with conditions: 

h(l) = 0, d y h(0)=0, A(0) = fi(0) = 0, / hdy = a- 1 . (101) 



J o 

Considering the evolution equation in (|41ap as y — > 1, we find a = B(l) as expected, where we 
have used the requirement that (hdyK} (1) = to alleviate the logarithmic pressure singular- 
ity, and that we also do not want a singularity in A at the contact line. This confirms Ucl = cl. 
To further simplify the equations for our numerical scheme, we note that we may rewrite A 
and B to avoid the calculation of the fifth derivative of the free surface. Consequently, we 
make the substitutions 

A = A — — P y —, B = /3B — - y , 102 

a + 12(1 + 4a)hp a + 12(1 + 4a)hp 
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Figure 11: The effect on cq when varying Rof3f for selected values of the other parameters 
using the finite contact line pressure condition. 



where a = 1 + 24a + 16a 2 for simplicity, so that the governing equations become 



dh ay dh d 
dt a dy dy 



ZAh 2 Bh (3/3 2 /i 2 + 4a)(l + 4a)/i 2 + a/3/i 3 d 3 /i 
2a + Ja~ + [12(1 + 4a)fi 2 h + a^]a 4 dy 3 



^ _ fpad^A + fp(l + 4a) d 2 B 



6/3/i + 1 + 4a h d 3 h 



(3a 2 dy 2 



4/3a 2 dy 2 a + 12(1 + 4a)h(3 a 3 dy 3 ' 



(103a) 
(103b) 



A 



T 

Jo 2 



dh OA 12/i/3 + 1 + 4a d 2 A d 2 B 



^ dy dy ^ 



dy 2 



+ 



■ d 2 h 



, + 3f3A-—r 
dy A dy z 

2 3(1 + 4a)(3h + 8a d 3 h 



with at y = 1 



B 



1 + 4a - Spa 
2p(l - 4a) 



+ 



A, 



a 3 a + 12(1 + 4a)h(3 dy 



h—r, (103c) 



dh\ 



ar d r 



dyj 



2/3 dy 



[1 + 4a - Spa + 12(3h) A + 2(p(l + 4a) + 2)B 



(104) 



and 



(/i<9 3 /i) (1) = 0, /»(1) = 0, d y h(0) = 0, A(0) = 5(0) = 0, / hdy = a~ 1 , B(l) = a/3. 

Jo 

(105) 

For a given free surface profile at a particular point in time, we may solve for A and B 
using the second two governing equations (|103bj) -( |103cp and boundary conditions (|104p . along 
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Figure 12: Velocity plots in the region close to the contact line of the two formulations of 
the interface formation model, (a) Xm = and (b) Xm = 1, in the frame of reference of 
the contact line. Parameter values chosen as (3 = 10 3 , a = 1/12, p = 0.1, (3f = 1, with 
additionally = 1/3 in (b). Plots shown in outer variables for perceptibility and to relate 
to Fig. [TOj and both use the finite contact line pressure condition. 



with A(0) = B(0) = 0. We then have the evolution equation (|103a|) . which contains a 
fourth order spatial derivative of h as in the Navier-slip case, along with the remaining 
five conditions of (|105p . Comparing the conditions with the Navier-slip case, we have here 
{hdyK) (1) = and a = B(l)//3, whereas for Navier-slip we would have (d y h) 2 (1) = a 2 , and 
a = (3nscl~ 3 (hdyti) (1). This would suggest that the additional condition, not listed above, 
is not required when enforcing finite pressure at the contact line — in agreement with our 
asymptotic analysis. 

To guarantee mass conservation we employ the transformation of Savva and Kalliadasis^i 
to solve for the integral of the free surface. We consider 



G = [ V h(y')dy', (106) 



so that the PDE (|103a]l becomes 



dG _ a 
dt a 



dG 
dy 



3Ah 2 Bh (3f3 2 h 2 + 4a)(l + 4a)/j 2 + a/3/i 3 d 3 h 



(107) 



2a [12(1 + 4a) /Ph + afla 4 dy 3 ' 

where we have integrated the evolution equation with respect to y. Now at y = 1 

G(l) = l/o, d y G{\) = 0, B(l) = Pa, (108) 

where the value of G(l) comes from the definition of G in ([1060 and the volume condition in 
(jlOip and we then automatically satisfy 

f My = G(l) - G(0) = l/o. (109) 
Jo 

We should point out a difference in the setup of the problem in comparison to the Navier-slip 
model. The contact line velocity condition a = B(l)//3 may still be applied here, whereas 
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the equivalent condition in Navier-slip is of the same order as in the PDE, and can no longer 
be used. Instead, we lose (/i<9 3 /i) (1) = in a similar manner, and then have the correct 
number of boundary conditions. The numerical solution of the evolution equation is based on 
spectral differentiation in space and adaptive, semi-implicit time stepping, following similar 
ideas from the scheme outlined in the Appendix of the study by Savva and Kalliadasis.— We 
note that along with the droplet radius behaviour we are also interested here in the variation 
of the microscopic contact angle, which given its dependence on the contact line velocity, is 
very sensitive to the initial condition imposed. Here, we specify an arbitrary droplet profile 
which has agreement with the leading order outer solution and with boundary layers of width 
0(ft~ 1 ), which is then allowed to briefly relax towards a quasistatic regime to provide the 
initial condition imposed. This then allows for a fair comparison between asymptotic and 
numerical results by minimising the impact of the specified initial condition. 

Figures [TBHTUI show the evolution of the droplet fronts and the microscopic dynamic con- 
tact angle behaviour in some instances for the Xm = interface formation model. Solid lines 
are from the full numerical solution of the partial differential equations, with dashed lines 
based upon the asymptotic results of the spreading rate in ([68|) . with cq determined through 
the solution of the ODEs in Sec. l3.2.2l for the droplet radius, and with the microscopic contact 
angle a~ 1 d y h\ y= i = —(1 + Pq). Unlike for the Navier-slip case,— when considering spreading 
outside the the region of asymptotic validity (given in (|7ip ). we are not always able to ob- 
tain agreement between asymptotic and full numerical solutions. As such, the initial droplet 
radius is always chosen such that (|71 j) holds. 

Figure [131 has parameters chosen to be /3 = 10 3 , ftf = 1, a = 1/12, p = 1 as in Figs. [5HS] 
and the streamline plot in Fig. [9]Ja) to show the excellent agreement between the asymptotic 
and full numerical results, where these typical parameter values have cq = —0.1 and Pq/cl = 
0.54 from our asymptotic results. Figure HH has the extreme parameter values ft = 10 3 , 
ftf = 100, a = 1/12, p = 10, which give cq = 56.96 and Pq/cl = 55.4. Excellent agreement 
is seen for the droplet radius, and good agreement for the microscopic contact angle for 
larger times. We would not expect the microscopic contact angle comparisons to show perfect 
agreement, especially for these extreme parameter values, due to its dependence on the contact 
line velocity and thus the initial condition prescribed. This can also be observed with slip 
and precursor film models, see Fig. 1(c) of Ref. [ljj 

Figure [15] gives two further examples of the droplet radius evolution using the parameters 
of the streamline and velocity plots in Figs. E^b) and 112(a). being (a) a(0) = 1.3, ft = 10 3 , 
ftf = 0.1, a = 1/12, p = 10 and (b) a(0) = 1.2, ft = 10 3 , ftf = 1, a = 1/12, p = 0.1. 
Finally, Fig. [TU1 demonstrates a situation for a receding droplet. The initial droplet radius is 
a(0)=2.3, with parameter values chosen as ft = 10 3 , ftf = 1, a = 1/12, p = 1, corresponding 
to c = -0.1 and P /a = 0.54. 



5.2 Numerical solution of the full problem for \ m = 1 

We now return to the governing equations of the full boundary value problem for x-m = 1> (|95p 
and (|34p . (|36p . (|38|) - ([39p . Transforming the system to the non- moving coordinate y through 
x = ay, the governing equations become 

3h ay 8h , 1 d ( h 3 d 3 h 3 A , 2 



m .^iS W + 2* +M|=w (110) 
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(a) 1.8 



(b) -0.99 




Figure 13: Interface formation model with Xm = 0: Evolution of (a) the droplet fronts, and 
(b) the microscopic contact angle for initial droplet radii a(0) = 1.4, with parameter values 
are chosen to be j3 = 10 3 , /3f = 1, a = 1/12, p = 1. Dashed lines use the asymptotic equations 
from matching with cq = —0.1 and Pq/cl = 0.54. Solid lines give the solution of the PDEs, 
and in (a) the lines are nearly indistinguishable. 

with the remaining equations as for Xm = 0, given in (|41bp - (|43p . and with the modified contact 
line velocity Ucl = B — aRo (c^/i) -1 D. Considering (jllOp as y — >• 1, and implementing the 
second boundary condition, we require 



which confirms Ucl = o,- The procedure follows that of the Xm = case closely, in particular 
the substitutions ()102j) and (|106p . and the conditions (jlOip . are employed, noting that ([3T]) 
and implementation of the contact line conditions gives the modified area condition: 



The mass of the droplet remains constant, but the area varies in the above condition as 
the proportion of the mass in the surface layers varies due to the mass transfer terms in 
the governing equations. Figures [T7] and [18] show excellent agreement between asymptotic 
and full numerical solutions for the evolution of the droplet fronts, with the microscopic 
dynamic contact angle behaviours also shown. Initial droplet radii are chosen to ensure the 
conditions of our asymptotic analysis are satisfied. Figure [T71 depicts the evolution using the 
parameter values of the velocity plot in Fig. [T2Tb). being f3 = 10 3 , a = 1/12, p = 0.1, f3f = 1, 
Ro = 1/3, where Co = 0.12 and Pq/cl = 0.093. Figure [18] uses parameter values /3 = 10 3 , 
/3f = 10, a = 0.1, p = 10, Ro = 0.1 to compare results for the larger value, Co = 10.44 and 
P /a = 12.09. 

6 Conclusions 

We have conducted a critical analysis of the theory of contact line motion proposed by 
Shikhmurzaev, in both original— and modern— formulations through application to the prob- 
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Figure 14: Interface formation model with x-m = 0: Evolution of (a) the droplet fronts, and 
(b) the microscopic contact angle for initial droplet radius a(0) = 1.65, with parameter values 
chosen to be ^ = 10 3 , /3f = 100, a = 1/12, p = 10. Dashed lines use the asymptotic equations 
from matching with cq = 56.96 and Pq/cl = 55.4. Solid lines give the solution of the PDEs. 



lem of a thin two-dimensional droplet spreading on a flat solid substrate. We confirm a 
number of features of the model in this scenario, in particular: 

• The interface formation model is able to alleviate the logarithmic pressure singularity 
at the contact line, a problem encountered with the Navier-slip model. (The Rucken- 
stein and Dunn slip model displays nonsingular pressure behaviour, as do other models 
mentioned in the Introduction, Sec. [TJ such as those with diffuse interfaces or precursor 
films). 

• For this particular flow scenario the model needs no adjustment. The requirement 
of finite pressure at the contact line closes the system and no additional contact line 
conditions are required. 

• If the logarithmic pressure singularity is allowed (given that it will still represent a 
finite force), then one may still adhere to the additional condition of Bedeaux— and 
Billingham,— but with a further model parameter to be imposed, alongside the myriad 
of parameters already necessary in the interface formation model. 

• The microscopic contact angle is determined through the solution of the governing 
equations, with no assumptions required on its velocity dependence from experimental 
results. 

It is instructive to make some further remarks about the model, relevant to both situations 
of enforcing finite pressure at the contact line or instead applying the additional condition: 

• For quasistatic droplet spreading in the long-wave approximation the interface formation 
model effectively reduces to a sophisticated slip model, with spreading rates based on a 
parameter Co, obtained through our matching procedure. 

• An analytical formula for cq in terms of the model parameters has not been possible, 
requiring instead a numerical solution of a system of two ODEs. 
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Figure 15: Interface formation model with Xm = 0: Evolution of the droplet fronts for initial 
droplet radii and parameter values (a) a(0) = 1.4, j3 = 10 3 , f3f = 0.1, a = 1/12, p = 10 
and (b) a(0) = 1.4, fi = 10 3 , j3f = 1, a = 1/12, p = 0.1. Dashed lines use the asymptotic 
equations from matching with (a) cq = 1.28 and (b) cq = 0.16. Solid lines are from the 
solution of the PDEs and are nearly indistinguishable. 

• There is a large number of parameters available to vary when compared to slip models. 
Even in the simple setting considered here, where assumptions have been made to keep 
them to a minimum (for instance that the gas-liquid and solid-liquid interfaces share 
the same values for r, 7, a and (3). Uncertainty in their values, especially if the model 
is to be applied to a range of different fluids in a variety of technological process, makes 
predictions difficult. Detailed and systematic measurements of the parameters, either 
experimentally or through molecular dynamics simulations, are needed. 

• The complexity of the model, even when applied to the simple situation adopted here, 
may be part of the reason why there is not a larger body of research using it — along 
with the questions about the physical basis of the model, discussed in Sec. [TJ As many 
of the crucial features of the model have remained in the long-wave approximation 
for the droplet spreading problem, the latter could serve as a useful starting point to 
consider further details of the model. On the other hand, the Ruckenstein and Dunn slip 
model seems to capture all the necessary physical ingredients of contact line dynamics 
with the exception of nanoscale rolling motion (the experimental verification, or indeed 
prohibition, of which is beyond the capability of current results known to the authors). 
It could thus be used as the basis for the development of models which capture this 
effect also, if desired, and yet are simpler and with a smaller number of parameters to 
the interface formation model. 

Taking full account of the mass transfer between bulk and surface layers, as in the modern 
formulation of the interface formation model, is required for nanoscale rolling motion through 
the contact line. The original formulation predicts a low-velocity region as seen for slip 
models, which prevents the fluid reaching the contact line in finite time. The evolution of the 
droplet radius using the interface formation model reduces to an equivalent expression for a 
slip model, with the microscopic dynamic contact angle having a correction to its static value 
proportional to the contact line velocity. This is also equivalent to the spreading rate for a 
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Figure 16: Interface formation model with Xm = 0: Evolution of (a) the droplet fronts and (b) 
the microscopic contact angle for initial droplet radius a(0) = 2.3 to demonstrate a receding 
droplet situation with parameter values j3 = 10 3 , ftf = 1, a = 1/12, p = 1. Dashed lines use 
the asymptotic equations from matching with cq = —0.1 and Pq/cl = 0.54, with solid lines 
from the solution of the PDEs. 



slip model with microscopic dynamic contact angle equal to the static contact angle, but with 
a slip length modified by a parameter cq, appearing in the outer expansion of the inner region 
equations. 

Our study has been based on the simplest of settings, that of a droplet spreading on a 
planar horizontal substrate, where the dynamic behaviour of the interface formation model 
could be investigated, both analytically and numerically. The addition of substrate topog- 
raphy, chemical heterogeneities, gravitational effects and the extension to three dimensional 
droplets to allow for comparison with experiments would all be of interest (albeit increasing 
the complexity of the analysis), as such effects have been shown to induce many interesting 
and often surprising phenomena (see e.g. the recent studies in Refs. Hfl. The value 



of the relaxation time t has been one of the points of contention for the interface formation 
model.— Provided an extension to three dimensional droplets produces similar results, an 
experimental comparison could be drawn for the evolution of the droplet radius. For known 
values of the other parameters, then cq, the constant determining the rate of spread of the 
droplet, could be fitted to the experimental data which in turn could be used to determine 
t. It would then be of interest to compare this to the range of values suggested by Blake and 
Shikhmurzaev,— and motivated the typical values adopted here and by Billingham. 24 i 37 
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Figure 17: Interface formation model with x-m = 1 ; Evolution of (a) the droplet fronts, and 
(b) the microscopic contact angle for initial droplet radius a(0) = 1.4 and parameter values 
(3 = 10 3 , /3f = 1, a = 1/12, p = 0.1, Ro = 1/3. Dashed lines use the asymptotic equations 
from matching with cq = 0.12 and Pq/cl = 0.093, with solid lines from the full numerical 
solution. 



A Simplified dimensional equations 



We record the simplified equations relevent to Sec. 12.11 having eliminated excess variables. 
On the solid-liquid interface we have 

7 dps 



dU dV 
~dZ + ~dX~ 
dp s s d 
~df~ + dX 



2p dX 
'p s s U\z=o 



P 



V = X r 



jap s 



dp; 



dX 



PS ~ PSe 

I ! 

PT 

PSe ~ PS 



(113) 



having eliminated u s s and ag using the last two equations of (|17j) . On the gas-liquid interface 
we eliminate (Jq using (20g). It is then possible to eliminate both components of the surface 
velocity u s G and v s G using (|20fjl and (|20dp to find 



U-j 



1 + 4a dp s G p G - p Ge dH 



4f3n dX 



X. 



. XT l + 4adp s c dH Pr ~ Pee 



4/3n dX dX 



then giving the surface equations 

dT +u dx~ v +Xm t P 
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P G - P + n G • T • n G = 7 

+ 4 

CA J 
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:au dv 
&- 2 ^dz + dx 



P(0)G - Pg d 2 H 

dU dH _ -/Vndp G 
~3X~dX ~ p ~dX' 



rpVU ' 
(114) 
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(115b) 
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Figure 18: Interface formation model with x-m = 1 ; Evolution of (a) the droplet fronts, and 
(b) the microscopic contact angle for initial droplet radius a(0) = 1.6 and parameter values 
/3 = 10 3 , j3f = 10, a = 0.1, p = 10, Rq = 0.1. Dashed lines use the asymptotic equations from 
matching with cq = 10.44 and Pq/cl = 12.09, with solid lines from the full numerical solution. 



where • T • is as given in (|19p . Finally at the contact line we have 



u - 



7 (1 + 4a) dp s G 
4/3% dX 



v PG - PGe dH rr 
Xm 77^7 ~ U CL 
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1 rpVn 



sm( 
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cos 9 d + p s G 
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V 



7(1 + 4a) dp s G dH 



(P(0)G ~ Pa) cos °d~ PS= (P(o)G - PGe) COS fl s - p| e , 

U dp s s 
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7a 



^o(p|-P5e + PG-PGe)> 



(116a) 
(116b) 
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having used a s Ge = 7(P( )G — Phe) ano - ^ ne equivalent in the solid surface. This completes our 
description of the governing equations for droplet motion, being able to eliminate the surface 
velocities and tensions, leaving a simplified set of equations. 
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